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Abstract 


The paper presents a new computational model of non-steady operation of a PEM fuel cell. The model is based on the macroscopic hydrodynamic 
approach and assumptions of low humidity operation and one-dimensionality of transport processes. Its novelty and advantage in comparison with 
similar existing models is that it takes into account the finite-time equilibration between vapor and membrane-phase liquid water within the catalyst 
layers. The phenomenon is described using an additional parameter with the physical meaning of the typical reciprocal time of the equilibration. A 
computational parametric study is conducted to identify the effect of the finite-time equilibration on steady-state and transient operation of a PEM 


fuel cell. 
© 2007 Elsevier B.V. All rights reserved. 
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1. Introduction 


The polymer electrolyte membrane (PEM) fuel cells are 
promising power sources with numerous potential and already 
realized applications. Although a subject of scientific scrutiny 
for several decades, the operation of these devices is still not fully 
understood. The main reason is the complexity of the underlying 
physical phenomena, which, in particular, concerns the intercon- 
nected processes of transport of gaseous and liquid species and 
phase transformation of water. 

It has been long understood that the experimental approach 
alone is insufficient and should be complemented by theoret- 
ical and computational models, which provide needed insight 
into the cell operation and possibility of the cell optimization. 
Numerous models have been developed over the last two decades 
(see Refs. [1,2] for a recent review and Refs. [3—5] for fur- 
ther examples). The past modeling studies identified the effects, 
which are vitally important for the fuel cell operation and require 
accurate modeling. A particularly important area, which can 
be greatly assisted with more accurate modeling, is the water 
management of fuel cell. The goal of the water management 
is two-fold. The membrane hydration has to be maintained at 
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a desired level. On the other hand, condensation of excessive 
water into droplets within gas diffusion and catalyst layers has 
to be avoided since it leads to blockage of the transport of gas 
reactants. 

In the present paper, we utilize the principles of macroscopic 
mass transfer to develop a model of a single fuel cell. The novelty 
of our approach is in the description of the transport and phase 
transformation of water within the catalyst layers. As opposite 
to almost all previous studies (the only exception [6] is dis- 
cussed below), we do not assume equilibrium between the pore 
vapor and membrane liquid water but rather consider them as 
different phases equilibrating at a finite rate. Let us introduce 
the phenomenon and discuss it in some details. 

Catalyst layers, the most complex and important parts of any 
fuel cell, generally consist of three phases: a solid matrix, a 
polymer membrane phase with attached agglomerates of liquid 
water, and pores filled with a gas mixture. The concentrations of 
liquid water retained in the agglomerates and water vapor in the 
pores are determined by their diffusive and convective transport 
and phase transition into each other. The traditional theoretical 
approach is to assume that the liquid and vapor phases within the 
catalyst layers are at equilibrium with each other described by the 
experimental uptake curve. The curve was originally obtained 
in Ref. [7] as an equilibrium relation between the average liquid 
water content within the membrane and the vapor concentration 
in the surroundings. It was mentioned in Ref. [7] that the equi- 
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librium was not achieved instantly but over a prolonged period 
of time determined by the slow diffusion of liquid water through 
the membrane. 

It has been long understood that the slow equilibration 
between the membrane water content and vapor concentration in 
anode and cathode has to be taken into account by any unsteady 
fuel cell model. This was done, for example, in the recent work 
[8]. The approach was, however, not followed in the description 
of the processes within the catalyst layers, where an equilibrium 
between vapor and liquid water was assumed. As it is discussed 
in Section 2.2 of the present paper, the assumption, although 
acceptable as a first approximation, is not entirely correct. 

The present paper questions the validity of the equilibrium 
assumption applied within the catalyst layers by considering 
that the ionomer water forms agglomerates and, thus, not all 
hydrophilic elements of the ionomer are in contact with the 
pore vapor. The complete equilibration would occur only over 
the relaxation time determined by the diffusion of liquid water 
through the agglomerates. Moreover, the combination of the 
finite-rate relaxation and diffusive water transport within the 
agglomerates and pores may lead to absence of equilibrium even 
in steady-state behavior. 

Addressing this issue, we incorporate the phase transition into 
the model of a single fuel cell and investigate the sensitivity of 
the cell performance to the absorption rate of vapor by ionomer 
in the catalyst layers. We use the approach of Ref. [6], where 
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a simple phenomenological model of the process was proposed 
and analyzed for the case of a separate cathode catalyst layer. 
The rate of equilibration is defined as a constant t, ! where Ty 
is the typical relaxation (equilibration) time, which primarily 
depends on the volume fraction of the ionomer in the catalyst 
layer and the area of the contact surface. The upper bound of 
this time scale can be estimated for the cell properties used in 


the present study (see Table 1) as T = (m" 3L /Dm.eff ~ 
0.5s. Here Le is the thickness of the catalyst layer, m is the 
volume fraction of ionomer in the catalyst layer, and Dm,eff is 
the effective diffusivity coefficient of the ionomer water related 
to the membrane diffusivity Dm via the Bruggeman correlation. 
The estimate should be understood as an unrealistic upper bound, 
since it is based on the assumption that a single agglomerate has 
the dimension L, of the thickness of the catalyst layer. The actual 
relaxation time ty is likely to be considerably, perhaps two or 
three orders of magnitude, smaller than ty". 

Since our primary interest is in the automotive applications, 
where the load and external conditions (temperature, air humid- 
ity, etc.) may vary with time, the focus of the investigation is on 
the transient regimes of the fuel cell operation. It is necessary to 
mention that, with few exceptions, the fuel cell models reported 
in the literature are limited to steady-state or quasi-steady-state 
solutions. 

The non-steady cell operation is governed by very different 
time scales. It is interesting to estimate them and compare with 


Table 1 

Physical, electrochemical, transport properties, and geometrical parameters of the fuel cell 

Quantity Unit Value 

Atmospheric pressure, Patm Pa 101,325 

Pressure of saturated vapor, Psat Pa log 9(Psat/Patm) = —2.1794 + 2.953 x 10-2(T — 
273.15) — 9.1837 x 1075(T — 273.15)? + 1.4454 x 
10-7(T — 273.15) 

Diffusivity of hydrogen in gas mixture, D¥2 m?s~! 1.1 x 1074(7/353)>/*(patm/P) 

Diffusivity of oxygen in gas mixture, D°? m?s~! 3.2 10~4(7/353)>/? (patm/P) 

Diffusivity of vapor in gas mixture, D42° m? s7! 7.35 x 1075(7/353)*/?(patm/P) 


Gas viscosity, u 


Open circuit potential, Vo 
Membrane diffusivity, Dm 


Tonic conductivity of the membrane, « 

Electro-osmotic drag coefficient, na 

Exchange current densities at the anode and cathode sides x reaction 
surface area, jê pand jep 

Reference hydrogen and oxygen, concentrations, C Ha and c% 

Transfer coefficients at the anode and cathode sides, a + @ and œe 

Dry membrane density, pary 

Equivalent weight of ionomer, EW 

Porosity of the anode, € 

Porosity of the catalyst layer, € 

Volume fraction of ionomer in catalyst layers, m 

Permeability of the anode and cathode, K 

Permeability of the catalyst layers, K 

Pressure at the inlet and outlet ends, po 

Thickness of anode and cathode, Lı and L5 

Thickness of the anodic and cathodic catalyst layers, L2 and L4 

Thickness of the membrane, L3 


kg (ms)! 9.88 x 1076C + 1.22 x 1075C™20 + 2.3 x 1075CO2 + 
2.01 x 1075C2 
v 1.23 — 9 x 1074(T — 298) 
m a 3.1 x 1077A (e0284 _ 1)e72346/T | A<3 
4.17 x 1078A(1 + 16le*)e~249/T, > 3 
A (Vm! (0.51394 — 0.326) exp[1268((1/303) — (1/T))] 
1.0 forà < 14 
Am? 10° and 104 
mol m~’ 40 
2 and 1 
kg m~? 1980 
kg mol! 1.1 
0.6 
0.4 
0.26 
m? 107!? 
m? 10715 
Pa 2Patm 
m 3 x 1074 
m 1075 
m 5.1 x 107-5 


All data are from Refs. [10,11]. 
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the relaxation time introduced above. The first is the typical 
diffusion time of the reactant gases through anode and cath- 
ode, which can be estimated as tg = L? / Dg ett, where Lg is the 
thickness of the anode or cathode gas diffusion layer and Dg eff 
is the effective diffusion coefficient of gaseous species. Typi- 
cally, tg = 0.001 — 0.05 s. The second time scale is related to 
the membrane absorption (desorption) of water. The process is 
characterized by a much larger diffusive time scale, which can 
be estimated as Tm = i /Dm, where Lm and Dp are the mem- 
brane thickness and diffusion coefficient, and is typically equal 
to 1-10s. 


2. Model 


A fuel cell consists of porous anode and cathode diffusive 
layers, membrane, and two catalyst layers of small but finite 
thickness. Considering the focus of the present work, the model 
does not include the anode and cathode gas channels, and the 
boundary conditions are applied directly at the outer borders of 
the anode and cathode. 

The mathematical model closely follows the model used in 
Ref. [2], with the major difference made in consideration of the 
catalyst layers. Low humidity regime of the cell operation is 
considered. The diffusion layers are assumed to have no liquid 
droplets. The vapor penetrating into the fuel cell at the anode 
and cathode sides is to be under partial pressure that is less than 
the saturated value, so no condensation occurs in these layers. 
The membrane (including the membrane phase of the catalyst) 
is filled with liquid water, and the phase transition between 
liquid water and vapor can occur at the surface of the mem- 
brane phase. Under the conditions of low humidity operation, 
the main task of water management is to prevent the membrane 
dehydration. 

The low humidity assumption does not impose any serious 
limits on the scope of this investigation. It was experimentally 
shown in Ref. [9] that the fuel cell can be successfully utilized 
at low humidification levels without using external humidifica- 
tion. Although smaller efficiency rates are attained under such 
conditions, there are also advantages of avoiding the use of the 
external humidifiers, such as reduction in the cell cost, mass, 
and volume. At present, the low humidity operation is consid- 
ered to be promising for small portable, including automotive, 
applications. 

For the sake of simplicity, a one-dimensional, transient 
approach is taken, with all fields being functions of the coor- 
dinate x across the cell and time t. The model can be easily 
generalized to two and three dimensions in future studies. Tem- 
perature of the cell is assumed to be constant and equal to 80°C. 

In the following subsections, we formulate the governing 
equations and boundary conditions separately for each part of 
the fuel cell. The behavior of the gas mixture within the porous 
media is governed by the mass balance for each species and the 
balance of total momentum. The distribution of the membrane 
water content is determined by the equation for the balance of 
mass. The values of the physical quantities used in the simula- 
tions are summarized in Table 1. 


2.1. Anode 


Anode is modeled as a porous medium filled with the gas 
mixture which contains hydrogen, water vapor, and sometimes 
nitrogen and other non-reacting gases. In the simulations pre- 
sented below, the transport of hydrogen and vapor is resolved 
directly, while the nitrogen concentration is found from the 
mass balance conditions. The mixture is transported through 
the porous medium by diffusive and convective mechanisms. 
The ideal gas model is adopted for the mixture. 

The mass balance is expressed by the continuity equation 


a. + | (ou) =0. (1) 
X 


Here € is the porosity of the matrix, u the filtration velocity and 
p the density of the gas mixture, which is defined through the 
molar concentrations and molar weights of components as 


p = MC + MOC™O 4 YyN2CN2, (2) 


The permeability coefficients of the electrodes are very small, 
so the conservation of momentum may be described by the Darcy 
law 


——u. (3) 


Here, p is the gas mixture pressure and K is the permeability of 
the porous media. 
Equations for species conservation are 


ə ə ə ə 

Cts. Cys (pyc ), 4 
“or + ax ) ox ( eft 3y os 
a H20 ð H20 a H20 a H20 

ad ea 20) — < | pþROL cO |. 
“ae + T ) ox Ta ©) 


The Bruggeman correlation is applied in order to find the effec- 
tive diffusion coefficients 


H 1.5 pH H20 1.5 pH20 
DR =e5pD®, DREO = e!5 DO, (6) 


Equation of state reads 
p = RT(C™ + C? + CN), ) 


where R is the universal gas constant and T is the constant 
temperature of the cell. 

At the outer (inlet) boundary of the anode layer, the concen- 
trations of hydrogen and vapor (humidity), i.e. C#2 and C¥2°, as 
well as pressure p are assumed to be known. The concentration 
of nitrogen can be determined from the equation of state (7). 

At the boundary separating the anode and anodic cata- 
lyst layer, continuity of hydrogen and vapor concentrations 
C™ and C#2° and continuity of their fluxes Dia ac? /ax, 
DIRO ac#29 /ðx are assumed. Continuity of pressure p and mass 
flux pu is also required. 


2.2. Anodic catalyst layer 


The material of the catalyst layer is a fusion of two fractions: 
a porous matrix filled with the gas mixture and the membrane, 
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where liquid water exists in the form of agglomerates. In the 
present work, both fractions are assumed to be homogeneously 
distributed over the thickness of the catalyst layer. The layer is 
considered as a porous medium with the pores formed by the 
solid matrix and water agglomerates. The porosity and perme- 
ability of the layer are constants, the values of which are reduced 
in comparison to the characteristics of the anode in order to 
account for the fraction occupied by the membrane phase. The 
transport within the catalyst layer is governed by the equations of 
the gas mixture filling the pores. Convection within the ionomer, 
as well as diffusion of gas components through liquid water are 
considered negligible. 

It is assumed that hydrogen oxidation reaction (H? — 2H* = 
2e7) occurs in the pores. Protons are then diffused though the 
membrane phase of the catalyst layer and the membrane to the 
cathodic catalyst layer. 

Traditionally, the membrane liquid water content (hydration) 
is characterized by the quantity à which is defined as the num- 
ber of water molecules per sulfonic-acid group (SO3_ ) (see, e.g., 
[7]). The liquid water content can be at equilibrium with the sur- 
rounding vapor at arbitrary vapor concentration. The situation is 
principally different from the simple evaporation—condensation 
process due to the additional internal energy of the membrane 
active centers that keeps the water molecules in the agglomer- 
ates. The equilibrium is defined by the uptake curve, traditionally 
expressed as a relation between the membrane hydration À and 
the vapor activity a = RTC H20 / Psat, Where psat is the saturated 
vapor pressure. According to Ref. [7], the uptake curve can be 
approximated as 


fr0<a<1 


forl <a< 3. 


(8) 


f 0.0043 + 17.81a — 39.85a* + 36.0a°, 
© ) 144+ 14(a— 1), 


Here the subscript e stands for the equilibrium. 

In earlier models that considered water transport within the 
catalyst layer, instant equilibrium (8) between the vapor and 
ionomer liquid was assumed. This assumption, albeit suitable 
as a first approximation, is not strictly correct. The membrane 
liquid water in the catalyst layer forms large agglomerates. Only 
a small portion of the ionomer adjacent to the agglomerate sur- 
faces is in direct contact with the pore gas. Any variation of the 
vapor activity is followed by immediate corresponding change 
of the liquid water content in the surface region but takes some 
finite time to diffuse into the interior. As a result, the equilibration 
is achieved only after certain typical time, which was estimated 
in Section 1 as 1073 to 107! s. Furthermore, there can be steady- 
states characterized by the absence of equilibrium between the 
vapor and liquid water. Finite rate equilibration is offset in such 
states by the spatial diffusion of water. 

It can also be mentioned that the instant equilibrium approach 
leads to an inconsistency in formulation of the boundary condi- 
tions at the catalyst borders. At the boundary between the anode 
and catalyst layer, vapor passes freely, whereas liquid water is 
trapped in the membrane phase and does not penetrate into the 
anode. At the opposite boundary, the assumption is made that the 
vapor does not diffuse into the membrane, while the liquid water 


of the catalyst membrane phase does. These boundary condi- 
tions cannot be modeled correctly if one-to-one correspondence 
between the vapor concentration and water content of membrane 
phase is assumed. 

The present analysis incorporates a more realistic descrip- 
tion of the water absorption and desorption by the membrane 
phase. The transition between the vapor and liquid water is 
described with the help of a phenomenological model pro- 
posed in Ref. [6] and conceptually analogous to the common 
evaporation—condensation model. The rate of change of concen- 
tration of liquid water and vapor is assumed to be proportional 
to the deviation of the actual hydration à from the equilibrium 
hydration àe given by (8). The phenomenological coefficient 
of proportionality is a constant equal to the typical reciprocal 
equilibration time t}. 

The governing equations in the gas pores include the conti- 
nuity equation for the gas mixture, Darcy law (conservation of 
momentum), and the diffusion equations for hydrogen and vapor 


dp ə ml, H20 H20 
= -—M™ M29y*(pmde — C°), (9 
E T g” oP! Vy" (PmAe 4 ) (9) 
a 
P Ëu, (10) 
ox K 
a a a 1 

H2 4 7) H2) — p% H2 ; 11 
ro ge ae ( eff 5 C ) 2F” ah) 
e2 cO 4 Ê cto) 
a a 

= Ê (pio? co) — (pare — CO), (12) 

ax V Ë ax A 


Here, F is the Faraday constant and j is the transfer current of 
the hydrogen oxidation reaction discussed below. 

The membrane phase is assumed to be impermeable to gases. 
The transport of liquid water is given by the diffusion equation 


ə a 50 9 AH nd , 
m~ Cp? = pe Ca ed 


+y*(Pmàe — OP). (13) 


The water transport by the electro-osmotic drag of water 
molecules by migrating protons is included in Eq. (13) with 
the assumption that electro-osmotic coefficient mq is constant 
[12]. 

Eq. (13), which did not appear in previous models, is required 
in our model since the equilibrium between the liquid water 
and vapor is not imposed. In addition to the terms describing 
the diffusive and convective transport and the consumption of 
hydrogen in the electrochemical reaction, the Eqs. (9)—(13) con- 
tain new terms describing the absorption or desorption of vapor 
by the membrane phase. C#2° = p,..a is the molar vapor con- 
centration, where a is the water activity. c0 = Pmà is the 
quantity characterizing the water content of membrane phase 
(Pm = Pary/EW, where pary is the density of the dry membrane 
and EW is the equivalent weight defined as the weight of the 
membrane per mole of sulfonic acid group). y* = Ty l is the 
proposed phenomenological parameter. In the following, we will 
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use the non-dimensional relaxation parameter y = tgy*, where 
Tg is the typical gas diffusion time, which is estimated as Ly D 
using the total thickness of the cell and the effective diffusivity 
of hydrogen at p = Patm. The parameter y is varied between 100 
(almost immediate equilibration) and 0.1 (very slow equilibra- 
tion). This corresponds to the range 11.4 < y* < 1.14 x 10fs7! 
and the dimensional relaxation time t, between 8.8 x 1075 and 
8.8 x 1077s. 

All diffusion coefficients in (9)-(13) are modified via the 
Bruggeman correlation as follows 


H _ -1.5 pH2 H20 _ -1.5 pH20 
Dat =e~D š Dat =€“ D ; 


DO, =m!” Dm, (14) 


where m is the volume fraction taken by the ionomer in the 
catalyst layer. 

The model also includes the equation of state and definition 
of the density of the gas mixture equivalent to (7) and (2). 

Additional equations are needed to describe the electric cur- 
rent within the catalyst layer. We assume that the membrane 
is impermeable to the electrons, while the solid matrix of the 
anode and cathode is ideally conductive. These standard assump- 
tions allow us to model the anode and cathode as equipotential 
surfaces 


®, = Oon the anode, ®, = Veen on the cathode. (15) 


Here, ®, is the potential of the electronic phase. In all of the 
following equations, ® stands for the potential of the ionic phase. 

The transfer current, j= ði/ðx, is given by the 
Butler—Volmer equation [2], which is used in its linearized form 


4 ‘ — a F Qe 
= j (C™12o, 1/2 Ma 
J= jal ) RT 


ja = Pap(CH) F. a6) 


The material properties jê p, c, Œa, and a, are given in Table 1. 
The distribution of potential within the catalyst layer is 


defined by the equation for charge conservation: 


0 a@ . 
ae (i) +j=0. (17) 
Here, Keff is the conductivity of the membrane phase of the cat- 
alyst layer, related to the membrane ionic conductivity « (see 
Table 1) via the Bruggeman correlation Keg = m!>x. 

The boundary between the anode and anodic catalyst layer 
was considered in Section 2.1. The additional boundary condi- 
tions for the fields c0 and @ are, respectively, that of zero 
flux of liquid water acO /dx = 0 and zero proton current 
d@/dx = 0. 

The boundary conditions at the boundary between the anodic 
catalyst layer and the membrane are based on the assumption that 
the gases do not permeate the membrane. This means zero flux 
of hydrogen and vapor dC? /dx = 0, aC#2°/ax = 0 and zero 
mixture velocity u = 0 (or, equivalently, dp/dx = 0). The liquid 
water content che and its flux D¥2G0.aC}20 /dx are assumed 
to be continuous. Furthermore, continuity is required for the 
potential ® and the proton current cd®/dx. 


2.3. Membrane 


The membrane is completely filled with liquid water (pres- 
ence of dissolved gases is neglected). The water transport is 
determined by diffusion and electro-osmotic drag. The latter 
mechanism is defined by the term proportional to —(dnqi)/(0x), 
where i is the current density and ng is the electro-osmotic coef- 
ficient. As reported in Ref. [12], na is close to 1.0 in the range of 
water content A < 14 considered in the present work. Since there 
is no chemical reactions within the membrane, i is constant there. 
The last two facts mean that the electro-osmotic drag, although 
present in the membrane, does not change the internal water dis- 
tribution. Owing to the osmotic drag, the water molecules are 
transported from the anodic catalyst layer to the cathodic one, 
which only affects the water content of the catalyst layers. 

The membrane is assumed to be almost impermeable for 
gases (the actual permeability is estimated as K ~ 107! to 
10720 mĉ), i.e., the filtration velocity u is assumed to be zero. 

The equation for diffusion of liquid water is 


0 ð 
mi (Ong cto) l (18) 
X 


The transport of protons in the membrane is described by the 
equation for charge conservation: 


ð op 


The boundary conditions at the anodic border of the mem- 
brane were discussed in Section 2.2. The boundary conditions 
at the cathode side are analogous. The only difference is that the 
flux of oxygen (instead of hydrogen) is assumed to be zero. 


2.4. Cathodic catalyst layer 


The gas species behavior in the cathodic catalyst layer is 
described by the model similar to that for the anode side. The dif- 
ference is in the composition of the gas mixture (oxygen, nitro- 
gen and water vapor) and in the type of chemical reaction, which 
is now the reduction of oxygen (2H20 — O2 — 4H* = 4e7). 

The equations of continuity, Darcy law, species conservation, 
and the liquid water transport are 


G) 1 
(ou) = -M — j 


MRO * x — Cibo , 
or ox 4F V (Pmàe — C3?) 
(20) 
dp u 
£ =-%u, 21 
ox K“ (a1) 
a a a a 1 
—C™ + —(uc™) = — | DB—c® , 22 
e a = ( eff gy 4F” (22) 
ð H20 ð H20 
EzE + ae ) 
ae pRO Ê co — y*(Pmàe — CO) (23) 
əx eff ax salts ~ « 
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ð 0 0 1/1 

— Cho ee p20. CH20 as 2 ots . 

ar = ( mef gg A raj 
+y*(Pmde — C129). (24) 


As in the anodic catalyst layer, a transport equation is used to 

describe the liquid water transport and source terms are added 

to account for the transition between liquid water and vapor. 
The equation of state is 


p = RT(C® + CŒ? + CN2), (25) 
The density of the gas mixture is defined as 
p = MC + MOCHO 4 YyN2CN2, (26) 


The distribution of potential within the catalyst layer is 
described by the equation of conservation of charge 


ð oP 


— — j= 0. 27 
TCAE (27) 


The transfer current in Eqs. (20)-(27) is given by the 
Butler—Volmer equation 


j = jeC™ exp(a®), (28) 


where the following dimensional parameters are introduced (see 
Table 1 for the physical properties) 

Ac F 
RT ` 


EENE: O2\—1 (5 ) _ 
Je = Sret(Cret) exp | [Vee — Vo) ), a = (29) 
RT 

The Bruggeman correlation (14) is used to find the effective 
diffusivity and ionic conductivity coefficients. 

The boundary conditions at the surface between the mem- 
brane and the cathode catalyst layer include zero fluxes of 
oxygen and water vapor 0C°2 /dx = 0, AC#2°/ax = 0, zero gas 
mixture velocity u = 0, and continuities of the liquid water 
content Cen, its flux D¥20 aCe? /dx, ionic potential ®, and 
current Kefp0P/ dx. 

The conditions imposed at the boundary between the catalyst 
layer and the cathode, are continuities of concentrations of gas 
species C02, C#2°, their fluxes D92 ac® /ðx, DER acH29 Jax, 
and the total mass flux pu. We also require zero fluxes of liquid 
water acho /dx = 0 and ionic current d@/dx = 0. 


2.5. Cathode 


The equations and the boundary conditions describing the 
mass balance within the cathode are identical to those for the 
anode and can be easily obtained by substitution of the quantities 
for the hydrogen by the corresponding quantities for oxygen. 


2.6. Method of solution 


The problem is simplified by eliminating the filtration veloc- 
ity u from consideration through substitution of the Darcy 
law into the continuity and species conservation equations and 
boundary conditions. The resulting system is solved by the 
finite-difference method. The spatial discretization is achieved 


by applying the second-order central discretization formulas to 
all spatial derivative terms written in the conservation form. The 
time integration is based on the explicit Runge-Kutta method of 
the second order. 

Each time step includes the actual time integration of the evo- 
lution equations and the solution of the boundary-value problem 
for the ionic potential. At the first substep, the time integration is 
performed, which results in new concentration fields satisfying 
the balance equations and boundary conditions. At the second 
sub-step, the updated species concentrations are used to find the 
distribution of potential at the new time level. For this purpose, 
the charge conservation equation is integrated over the catalyst 
layers and the membrane with no-flux boundary conditions at 
the outer boundaries. Solution of the boundary-value problem is 
determined by the shooting algorithm. Determining the potential 
field in this fashion is computationally intensive. The field, how- 
ever, varies significantly only when the membrane water content 
changes. One also has to take into account that the time step used 
in the present calculations is very small, several orders of mag- 
nitude smaller than the typical time of water diffusion through 
the membrane. It is, therefore, unnecessary to adjust the poten- 
tial at every time step. We found that, in the cases considered in 
the present paper, adjustment after every 1000 time steps is suf- 
ficient for accurate computation. The coefficients of diffusivity, 
conductivity, and viscosity are renewed at the beginning of each 
time-step. 

The finite-difference grid is uniform within each sub-domain 
of the cell (anode, cathode, membrane, and catalyst layers). The 
grid steps are different in different sub-domains. We conducted 
a series of computations with numerical grids consisting of 10, 
20, and 30 points in each subdomain. Numerical convergence 
was clearly observed for both transient and steady-state solu- 
tions. The solution fields obtained with 20 and 30 points in each 
subdomain were almost indistinguishable, which indicated the 
accuracy and acceptability of a 20-point grid. All the results pre- 
sented below in this paper were obtained using this grid. The time 
step was constant and equal to 4.4 x 107°? s, which constitutes 
5 x 1077 tg. 

As initial conditions, the uniform distributions of concentra- 
tion of hydrogen in the anode and of oxygen in the cathode, of 
vapor and nitrogen in both parts, as well as uniform pressure 
throughout the entire fuel cell were specified. The water content 
in the membrane was assumed to be initially in equilibrium with 
the vapor concentration. The cell was assumed to be fed by pure 
hydrogen at the anode side and by air (CN? /C92 = 0.79/0.21 
on the cathode side). 


3. Numerical results 
3.1. Model validation 


The model was validated by comparing the steady-state 
results with the results of three-dimensional model [11] cho- 
sen because of the closest agreement with our model in regard 
of the physical effects taken into account and availability of a 
detailed description of the simulation parameters. The assump- 
tion of the uptake curve equilibrium between the vapor and liquid 
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Fig. 1. (a) Comparison of polarization curves. Solid line represents the results obtained with the present model, the dashed line corresponds to the model [11]. (b) 
Stoichiometric coefficients defined by (30). External conditions correspond to 50% humidity of the fuel and air. y = 100. 


water in the catalyst layers was made in Ref. [11]. As a nearest 
approximation achievable in the framework of our model, the 
simulations used for comparison were performed in the limit of 
fast equilibration (large relaxation parameter y). We found that 
y = 100 was sufficiently large to represent the asymptotic limit 
y — œ, the solution being insensitive to further increase of y. 

Fig. la shows the comparison of the polarization curves. 
The curves are in good qualitative agreement, although there 
are some quantitative differences. Such a discrepancy can be 
explained as follows. First, there is an inherent difference 
between 3D and 1D models. It was shown in Ref. [11] that 
the multi-dimensional effects influence the cell’s performance. 
The distributions of concentrations and other fields vary signif- 
icantly along the cell. Any 1D model can only be considered as 
an approximation representing the bulk properties. 

Further explanation is related to the difference in the spec- 
ification of the external boundary conditions. In Ref. [11] and 
many other models, the boundary conditions are imposed at the 
inlets and exits of the gas channels. The gas species distributions 
were found to vary significantly along the channels. In our model 
the channels are not included and the boundary conditions are 
imposed at the outer surfaces of the anode and cathode diffuser 
layers. This does not allow us to reproduce the conditions of 
earlier simulations exactly. For example, the curve correspond- 
ing to the results of Ref. [11] in Fig. la was obtained at 50% 
and 0% humidity at the inlets to the anode and cathode channels, 
respectively. Actual humidity at the channel/diffuser boundaries 
varied significantly. Figs. 5 and 7 of Ref. [11] show the humidity 
in the anode channel varying between approximately 25% and 
75% and humidity in the cathode channel growing from 0% to 
about 100%. The distributions corresponding to different points 
of the polarization curve also differ from each other. In order to 
imitate the conditions of Ref. [11], we assumed 50% humidity 
at both the anode and cathode outer boundaries. 

Another aspect also related to the specification of the bound- 
ary conditions is that, in Ref. [11] and some other models, the 
inlet velocities of each species were specified by the values of 
the stoichiometric coefficients, defined as 
H2 O2 
Z 2C uE pa AC; ee (30) 


l l 


Ea 


Here ua, uc are the inlet velocities (the sign is taken so that 
positive ua Or Uuc correspond to a flow into the cell) and C tz ,C O2 
are the hydrogen and oxygen concentrations at the anode and 
cathode side, respectively. These coefficients define the ratios 
between the amounts of each reactant coming into the cell to 
the amounts of reactant consumed inside. These coefficients are 
always assumed to be more than 1 (e.g., Ea = 1.2 or 2 and & = 2 
in Ref. [11]). The present model, as mentioned earlier, does not 
include the gas channels. The model assumes that the flow of 
the reactant gases into the anode and cathode layers is induced 
by their consumption inside and, so, the entire amounts entering 
the cell are burnt inside. The values of pressure and species 
concentrations are specified at the surfaces of the diffuser layers 
as input parameters. The values are taken to be the same at 
any voltage and current density and, thus, are different from the 
current-dependent values that would be obtained in simulations 
including the gas channels and using the boundary conditions 
on the stoichiometric coefficients (30). 

As an illustration of our approach to formulation of bound- 
ary conditions, the stoichiometric coefficients estimated on the 
basis of the obtained solution are shown in Fig. 1b. The figure 
depicts that the coefficients are only weakly dependent on the 
cell voltage or current densities, which means the values of the 
inlet velocities vary proportionally to the values of the current 
density. It is worthwhile to point out that € and & represent 
only convective transport into the diffusive layers. The analo- 
gous coefficients representing the full (convective plus diffusive) 
transport are always 1 in our model. The data in Fig. 1b, thus, 
provide non-trivial information on the relative significance of 
the two transport mechanisms. a exceeds 1, which means that 
the diffusive transport is negative (out of the cell). On the cath- 
ode side, & is negative corresponding to the convective transport 
of oxygen out of the cell. The net positive transport into the cell 
is provided by the strong diffusive flux caused by the negative 
gradient of the oxygen concentration. 


3.2. Steady-state regimes 


The steady-state solution obtained for the point Veeg = 
0.65 V of the polarization curve is illustrated in Fig. 2. The 
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Fig. 2. Steady state solutions at Veen = 0.65 V, 50% humidity of the fuel and air, y = 100. (a) Vapor concentration in anode, cathode, and catalyst layers; (b—d) 
liquid water content, ionic potential (in units of open circuit potential Vo), and current density in catalyst layers and membrane. The membrane domain is indicated 
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distributions of vapor concentration C H20 membrane water con- 
tent A, ionic potential ®, and current density i are shown. The 
profiles of gas species concentrations and pressure (not shown) 
demonstrate very small variations (less than 5%) within the cell. 

Many features of the distributions shown in Fig. 2 are typical 
for 1D simulations. The humidity profiles in Fig. 2a form as a 
result of the balance between the electrochemical production of 
water within the cathode catalyst layer, osmotic drag, and the 
diffusive and convective transport in the presence of constant 
humidity conditions at the outer surfaces of anode and cathode. 
The water content A within the membrane shows almost lin- 
ear growth from the anode to cathode side. Such a distribution, 
albeit similar to the results of earlier one-dimensional models, 
is in a disagreement with certain experiments, such as [13]. It 
has to be stressed that the linearity was not presumed by our 
model. The distribution of à was obtained as a solution of the 
transport Eq. (18) with a A-dependent diffusivity coefficient Dm 
(see Table 1). The experimental distribution could be affected by 
other factors not taken into account by the model. The most sig- 
nificant among them are the possible variability of the osmotic 
drag coefficient ng [13], two- and three-dimensionality effects 
[5], and the pressure gradient [3]. We note that, despite possible 
inaccuracies in representing the space distributions, our model 
provides the volume-averaged values of 4, which are very close 
to the results of earlier works, such as 3D computations of Ref. 
[11]. 


The ionic potential and current density are shown in Fig. 2c 
and d. The rate of ohmic voltage drop is almost constant within 
the membrane. The change of ionic conductivity « due to the 
variation of liquid content À is insignificantly small to cause any 
appreciable effect on the ®-profile. Within the catalyst layer, 
Kept (= m!°«) decreases sharply, which results in increased 
ohmic losses in the regions near the membrane. 

A novel feature of the present model is that vapor and liquid 
water are considered separately from each other in the catalyst 
layers and are not required to be in equilibrium with each other. 
This opens an interesting question of the effect of the relaxation 
parameter y on the steady-state behavior of the cell. To answer 
it, we conducted simulations at Vee = 0.7 V and with y varying 
between 0.1 and 100. The results are illustrated in Fig. 3, which 
shows that the entire distributions of vapor concentration C#2° 
and liquid water content à change significantly with y. Com- 
putations performed for other points of the polarization curve 
showed that the effect becomes even more pronounced at larger 
current density. 

Fig. 3a shows that the value of y affects the vapor concen- 
tration C#2° in the catalyst layers. A decrease of y (increase of 
the relaxation time t,,) leads to the growth of C H20 in the anode 
catalyst layer and a decrease in the cathode catalyst layer. The 
overall effect is that of reduced gradient of vapor concentration 
across the anode and cathode. The impact of changing y on the 
liquid water content is also significant (see Fig. 3b). A decrease 
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of y results in a stronger variation of A across the membrane. 
In the catalyst layers, the distributions of à form clearly visible 
zones of boundary-layer behavior near the catalyst-membrane 
surfaces. Widths of the zones and the amplitudes of variation 
of A within them grow with decreasing y. The zones can be 
detected in our model since the assumption of the uptake curve 
equilibrium between C#2° and A is not made. It can be deduced 
that using this assumption would inevitably lead to non-physical 
behavior at the membrane-catalyst border, either violation of the 
no-flux condition for C#2° or discontinuity of A. 

Fig. 3c and d show comparison between the distribution of 
à across the catalyst layers and the equilibrium distributions Ae 
calculated on the basis of C#2° using the uptake curve approx- 
imation (8). The results show that the difference between the 
actual and equilibrium water contents never fully disappears 
even though A becomes closer to Ae with increasing y. We would 
like to stress here that the curves for y = 100 essentially repre- 
sent the limit at y — oo since a further increase of y does not 
lead to any discernible variation of the solution. 

These results also elucidate the effect of the relaxation param- 
eter on the polarization curve. We found that, at given Veel, 
the cell current is virtually insensitive to y. A likely explana- 
tion can be seen in Fig. 3b. The modification of the profile of 
A, albeit significant, does not lead to substantial change of the 
total water content, average hydration and, thus, the membrane 
resistivity. 


3.3. Transient regimes 


The dynamic response of the cell to a sudden variation of the 
voltage Veen is illustrated in Fig. 4. The results are obtained by 
considering step changes in cell voltage from 0.7 to 0.6 V and 
from 0.6 to 0.7 V. The transient evolution is characterized by very 
different time scales. The surface overpotential adjusts immedi- 
ately resulting in an immediate change of current density (Fig. 4a 
and c). This causes corresponding changes in the concentrations 
and fluxes of reactants and other gases in the pores of the dif- 
fuser and catalyst layers, which occur on the very small diffusive 
time scale Tg (about 107°? to 107? s as estimated in Section 1). 
Our investigation is focused on the evolution of water and vapor 
distribution fields. This process is characterized by two larger 
time scales, that of the water diffusion in the membrane phase 
and the new time scale t, of the equilibration between the liquid 
and vapor within the catalyst layer. 

The evolution of the average water content of the mem- 
brane phase of the catalyst layers in Fig. 4b and d indicates 
the existence of two time periods. First, lasting about 0.5s 
is characterized by relatively fast change of (A). Its dura- 
tion increases with decreasing y. We can assume that during 
this period the evolution is dominated by the equilibration 
between the vapor and liquid water. The second period is 
much longer and characterized by slow changes of (A) both 
in the catalyst layers and the membrane. This evolution is 
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likely to be related to the slow water transport within the 
membrane. 

The variation of the average water content (A) in the catalyst 
layers is rather large especially in comparison with the small 
variation of the same quantity in the membrane. A further obser- 
vation can be made that, in agreement with the results obtained 
for steady-state solutions and illustrated in Fig. 3, the evolution 
and final asymptotic values of (A) in the catalyst layers are sig- 
nificantly affected by the relaxation parameter y. The difference 
between the curves for y = 0.1 and 100 reaches 10%. On the 
contrary, the behavior of average water content of the membrane 
is only weakly sensitive to y. This translates into the weak sen- 
sitivity of the cell current. One can see in Fig. 4a and c that the 
change of y can be responsible for a modification in the late 
stage of the current evolution controlled by the hydration of the 
membrane and accounting for few percent of the total current 
variation. 

In another simulation experiment, we analyzed the dynamic 
response of the cell to a change of the humidity conditions. The 
humidity at the cathode and anode inlet was simultaneously and 
suddenly changed from 50% to 100%. Such an experiment does 
not, strictly speaking, represent normal operating conditions, 
although may occur under special circumstances such as failure 
of the humidifier. It is also useful since it allows us to investigate 
the effect of the finite equilibration rate on an unsteady behavior 
that directly involves change of the membrane water content and 
vapor concentration of possibly lower amplitude. 


A comment is in order regarding the use of the low humidity 
approximation in the 100% inlet humidity case. It was discussed 
in Ref. [2] on the basis of earlier calculations that even at cur- 
rents as high as 1.5A x cm7? the level of saturation of liquid 
water in the cathode diffuser layer does not exceed 5%. At lower 
currents, such as those considered in the present paper, the satu- 
ration is even lower and unlikely to appreciably affect the oxygen 
transport. 

Fig. 5 shows the evolution of the current density i and the 
average liquid water content (A) in the membrane and the cat- 
alyst layers. The current density clearly follows the membrane 
water content in agreement with (A) being the main factor of 
the ionic resistivity. The evolution of the water content can be 
divided into three stages. During the first very short stage (at 
the time scale tg), the new humidity levels are diffused through 
the anode and cathode. The dominant phenomenon of the sec- 
ond stage is the phase transition from vapor to ionomer liquid 
water within the catalyst layers. As a result, the water content 
of both catalyst layers increases (decreases) sharply as shown 
in Fig. 5b and d and becomes substantially larger (smaller) than 
the water content of the membrane. This causes slow diffusion 
of water through the membrane, which characterizes the third 
stage of the evolution. The diffusion occurs on the largest time 
scale of the system. It can be seen in Fig. 5 that the asymptotic 
steady-state levels of (A) are only achieved after about 8s at 
y = 100 and y = 10 and after substantially longer periods at 
smaller y. 
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Second and third stages overlap, which does not allow us 
to clearly identify the duration of the second stage and hinders 
the evaluation of the effect of the equilibration parameter y. 
The effect is, however, more significant than the case of the 
transient behavior caused by the voltage change. Lower y delays 
the evolution of the membrane hydration and, thus, of the electric 
current. One can see in Fig. 5a and c that the variation of i with 
y, albeit small at the final asymptotic states, is comparably large 
(up to 25%) during the transient phase, which lasts few seconds. 


4. Concluding remarks 


This paper presents a one-dimensional isothermal model of a 
low-humidity non-steady operation of a PEM fuel cell. The main 
novelty and advantage of the model in comparison with simi- 
lar existing models is that it does not rely on the non-physical 
assumption of the uptake curve equilibrium between the pore 
vapor and ionomer water in the catalyst layers. Instead, the 
transition between the two phases is modeled as a finite-rate 
equilibration process. This introduces the typical equilibration 
time as the third governing time scale of the problem in addi- 
tion to the typical times of the pore gas diffusion and membrane 
water diffusion. 

A linear model is used to represent the equilibration process. 
The rate of equilibration is determined by the phenomeno- 
logical parameter y, which has the physical meaning of the 


non-dimensional reciprocal equilibration time t,. The equili- 
bration time is determined by the typical size of the ionomer 
water clusters in the catalyst layer and is likely to be of the 
order of 1074 to 107? s. Since no reliable experimental data are 
available for the size of the clusters, we chose to vary y ina 
range between 0.1 and 100 corresponding to the equilibration 
time 8.8 x 1075 < Ty < 8.8 x 107? s and conduct a parametric 
study of the effect of y on steady-state and transient behavior of 
the cell. 

The study shows that the equilibrium between the vapor and 
liquid water is never fully achieved within the catalyst layers. 
Even in the asymptotic limit of fast equilibration, which is rep- 
resented by the case y = 100 in this study, the liquid water 
distribution is noticeably distinct from the equilibrium distri- 
bution given by the uptake curve relation. This is true for both 
steady-state and transient situations. 

The steady-state solutions are found to be affected by the 
value of y. The effect on the distributions of liquid water con- 
tent within the catalyst layers and the membrane is particularly 
strong. However, this does not lead to comparably noticeable 
variation of the cell current. This may be explained by consider- 
ing that the average water content and, thus, the total resistance 
of the membrane do not change significantly. Another interesting 
feature of the steady-state solutions is the boundary layer-like 
behavior of the liquid water content near the surfaces separating 
the membrane and catalyst layers. 
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The study also investigates the transient response of a fuel 
cell to a sudden change of the voltage or external humidity con- 
ditions. It is found that the finite-rate equilibration between the 
catalyst vapor and liquid water can be a significant factor in 
understanding the internal behavior of the cell. Most importantly, 
it results in appearance of the additional stage of the response 
characterized by rapid change of (A) within the catalyst layers. 
This stage precedes and influences the change of the membrane 
hydration. Its duration is of the order of t,,, which implies sig- 
nificant impact of the equilibration parameter on the internal 
cell evolution. The fact that no appreciable variation of the cell 
current with y is found in our simulations does not exclude such 
a possibility in situations with different operational conditions 
and/or different set of physical features included in the model. 

To conclude, the parametric study presented in this paper 
shows the importance of the finite rate equilibration within the 
catalyst layers for the processes of water transport in PEM fuel 
cells. It is highly desirable to incorporate this phenomenon into 
the models of cell behavior in order to gain fundamental under- 
standing of transport processes within the cell. However, since 
the finite-time equilibration does not have any significant effect 
on cell current, the inclusion of this effect may not be neces- 
sary for zero-dimensional models designed to predict the overall 
performance of a cell. Further investigations are needed, espe- 
cially those directed toward providing an accurate estimate of 
the equilibration rate. 
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